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Abstract. Star clusters and dwarf galaxies gradually dissolve as they move in the potential 
of their host galaxy. Once their density falls below a certain critical density (which is com- 
parable with the background density of the galaxy) it is often assumed that their evolution 



> 



in 
o 



is completed. In fact the remnant of such a system forms a distribution of stars which are 



unbound to each other and which move on similar orbits in their host potential. With this 
motivation we study the evolution of an idealised unbound system and follow its expan- 
sion and dissolution in the tidal field of a model galaxy. Initially the stars are uniformly 
distributed (with a density below the critical density) within an ellipsoidal volume. The sys- 

03 ■ 

tern itself travels on a circular orbit within a galaxy modelled as an isothermal sphere. The 



initial velocities of the stars are chosen by assuming that they move on (three-dimensional) 
epicycles with guiding centre at the centre of the ellipsoid, though the usual epicyclic theory 
is altered to account for the self-gravity of the system. This is believed to be the first exact 
equilibrium model of a stellar system in a tidal field. Our main task is to study the stability 
of such configurations and the time-scale of their dissolution, as a function of the initial 
density and size of the ellipsoid. If the time of dissolution is measured by an increase of 
the half-mass radius of 50%, we find that systems of low density (~ 1% of the background 
density) and small radius (50 pc on an orbit of radius 10 kpc) can survive for about 20 
galactic orbits. For small systems we show that the lifetime is approximately proportional 
to the inverse square root of the density. 

Key words, galaxies: star clusters - methods: N-body simulations - Galaxy: kinematics 
and dynamics 



1. Introduction 



In recent years, considerable effort has been spent in examining the destructive processes which 
lead to the dissolution of star clusters and dwarf galaxi es, and the influence of these processes on 
the distribution function of these objects (e.g. iBaumgardt 1 1 1 998 ) and on the host galaxy. There 
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are two main areas of interest. First there is t he d issolutio n caused by inter nal processes such 



as evaporation (e.g. 



Giersz & Heggie 



1997 



Baumaardt & Makino 



is the merging o f these objects with the body of the galaxy (e.g. 



with each other (Fellhauer et al. 



2003). The second issue 



Velazauez & White 



1999) or 



2002). Within the context of the first problem, investigations 



usually stop when the object has dissolved according to some tidal stability criterion (Sec.|2l>. So 
far little attention has been paid to their subsequent fate, i.e. what happens to these objects after 
they become unbound entities. 

Retrograde motions are expecte d to pla y an important role in these unbound systems. It has 



been known for a long time ( Henon 



1969) that stars in retrograde motion about a star cluster (at 



least for clusters on circular galactic orbits) are stable (i.e. orbitally stable) even outside the tidal 
boundary; their orbits are best thought of as perturbed epicycles. As a star cluster loses mass and 
its tidal radius shrinks, we may expect that it is surrounded by some stars moving approximately 
on such orbits. 

Based on these considerations, the present study starts with a very special toy model in which 
all the stars move on retrograde orbits resembling epicycles with a common guiding centre. Their 
spatial distribution is a homogeneous ellipsoid, and the centre of the ellipsoid moves on a circular 
orbit within the potential of an isothermal sphere, which is taken as a model of the host galaxy. 
The simplicity of our setup permits some analytical investigation which gives useful insights 
into the results of numerical simulations. It is also believed to be the first example of an exact 
equilibrium of a stellar system in a tidal field. 

Having set up the sequence of models, our main interest turns to their stability. In particular, 
we aim to investigate the time-scales on which they evolve and finally disperse, and their spatial 
and kinematic properties as they do so. 

In the next section we briefly go through the theory of epicyclic motion with modifications 
to account for the self-gravity of the system itself. We also explain the setup of our toy model. 
Then we describe the results of our simulations, which cover a large parameter space. Finally we 
give our conclusions, and discuss their possible relevance for real systems, such as those found 
in and around our Milky Way. 



2. Theory & Setup 

The basic idea of our model is the following. Inside an ellipsoid of uniform density in a 
linear tidal field, orbits are modified epicycles, with a harmonic motion in z of amplitude 
zo, say (c.f Fig.Q). We consider epicycles which are centred at the centre of the ellipsoid. 
For a given epicyclic amplitude, we construct a distribution of zq so that the space density 
of this superposition of epicycles is independent of z (up to the edge of the ellipsoid). Then 
we construct a distribution of epicyclic amplitudes so that the space density is also uniform 
in x and y. Because we have arranged for the amplitude of the epicyclic x- and y- motions 
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to be in the same proportion as the corresponding axes of the ellipsoid, we can build up the 
entire ellipsoid in this way. 

In principle this gives rise to an equilibrium model. Our N-body simulations, however, 
do not maintain equilibrium because of relaxation, because the orbits are dynamically un- 
stable, and because the iV-body simulations use the exact tidal field and not a linear approx- 
imation. 




Fig. 1. A schematic representation of a uniform ellipsoid and one orbit. 



Let X, Y, Z be galactocentric coordinates, with Z perpendicular to the plane of motion of the 
system under study, and let R 2 = X 2 + Y 2 . Epicyclic theory is a linearised approximation in which 
stars move in a potential 



(D = 



1 3 2 <£ eff 



2 dR 2 



(R-D) 2 + 



1 <9 2 O eff 



2 dZ 2 



( Binnev & Tremaine 



(1) 



D.O 



1987), where D denotes the distance to the axis of symmetry of the galaxy 



and O e ff (R, Z) is the effective galactic potential. We set 

3 2 O eff 2 

D,0 

= v 2 , 



dR 2 

eff 



<9 2 <t>, 



dZ 2 



(2) 
(3) 



D.O 



where k is the epicyclic frequency. For the potential of a singular isothermal sphere these fre- 
quencies are given by 

k = V2Q, (4) 
v = Q, (5) 

where Q = Vq/D is the angular velocity of motion on a circular orbit of radius D around the host 
galaxy and Vo is the circular velocity. 
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We adopt initial conditions in which stars are distributed uniformly within an ellipsoid with 
semi-major axes a,^a,a, where a > is a free parameter, and £ > 1 will be determined from 
the shape of closed epicycles. Therefore we have to take the potential of a prolate ellipsoid of 
uniform density p into account. We use rotating, cluster-centered coordinates x, y, z with origin at 
the centre of the ellipsoid of stars. The coordinate axes are directed away from the galactic centre, 
in a direction opposite that of the circular motion round the galaxy, and orthogonal to the plane 
of m otion, respectively. The appropriate differential equations in the epicyclic approximation are 



then ( Chandrasekhar 



1942) 



x + 2Q.y + (k 2 - 4Q 2 )x = -B\x (6) 
y-2Q.x = -B 3 y (7) 

z +v 2 z = -Biz. (8) 

Here, B\ = 2nGpA\ and B 3 = 2nGpA 3 , wher e Aj and A3 are given by st andard formulae for the 



potential of a homogeneous prolate ellipsoid ( Bin nev & Tremaine 



1987, Table 2-1). 



This additional acceleration alters the frequencies of motion and the shape of closed epicy- 
cles. Eigensolutions of these equations give frequencies k' and V defined by 

K a = - ((k 2 + B x + B 3 ) + VO 2 + B X - B 3 ) 2 + 16Q 2 B 3 ) (9) 
v' 2 = v 2 +B, (10) 



and the epicyclic ratio 

k' 2 + k 2 -Bi _ 2Qk' 
2Q.k' ~ k' 2 - £3 



?= ^ - ~a 7T (ID 



Note incidentally that the two expressions for £ in ea. il 1\ are identical if eq.@ and the upper 
sign in eq.Q are used. 

If \k\ » B\, B 3 , then eq.Q impl ies that k' — k (if we choose the upper sign), and eg. I ll li t gives 
approximately the familiar axial ratio of epicyclic motion. In the general case, since we also use 
£ for the axial ratio of the ellipsoid, the right side of eg. I ll 1> depends on £ . Therefore we have 
solved for k' and £ using the above equations iteratively (taking the upper sign in eq.(|5})- 

A second solution of the differential equations (eqs.l|6j and (0) is obtained if one takes the 
lower sign for k' . It leads to an imaginary frequency ik = k" provided that 

p < ~r> (12) 

where po denotes the density corresponding to the isothermal background. (This is always ful- 
filled in our parameter range, though there is a series of bound systems (i.e. systems not satisfying 
eg.l !12» where all frequencies are real.) We can see that k <k k' if p is much smaller than the crit- 
ical density (cf. eg. (ll2» . and then we have: 



k 2 = — ((k 2 + B x + B 3 ) - V(* 2 + Bi - B3) 2 + 8k 2 B 3 ) 



2 

B 3 (13) 
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2Qk 



£.' - 

k 1 + Bi 

n 

(14) 



in which f is the ratio of the components of the corresponding eigenvector. 

Within the epicyclic approximation the complete solution is a linear combination of both 
solutions, i.e. 

x(f) = jtocosfVf + i/0 + Aexp(kt) (15) 
y(t) = £x sin(Af'f + iff) + Ag exp(kt) (16) 
z(t) = z cos(v't + cf>) (17) 

where Xq, Zo, ft, <fi and A are free parameters. (As with any unstable equilibrium when there is 
no dissipation, there is also an exponentially decaying solution, which we have ignored. It corre- 
sponds to the stable invariant manifold of the epicyclic motion, and is of no practical importance, 
being of measure zero.) 

As already mentioned the starting configuration is a homogeneous distribution of particles 
in an ellipsoid with axes a, £a and a. To achieve the correct spatial distribution we choose the 
parameters xq and zo according to the following distributions 



F(x ) = 1 . 



(18) 



F(zo) = (19) 
a x o 

with xo € [0, a] and zo e [0, yja 2 - x^]. The parameters ft and <f> are chosen randomly in the 
interval [0, 2jt]. 

Our recipe for constructing this model is given by eas.dT5l-fP7l (with three similar equa- 
tions for the initial velocities and adopting A - for the initial setup) and eas.dl8l and fl9l . 
Here we show that this gives rise to the desired uniform space density. 

From eqs.QD and fi"9t it follows that the joint probability density of the positive quan- 
tities xq,zo is 

ei \ 3xoZ(), 2 2 2\-l/2 /on\ 

f(xo,zo) = —^- (a - x o~ z o) ' ( 2 °) 

over the domain x > 0, zo > 0, x 2 + z\ < a 2 . We also take the distributions of ft, tf> to be 
independent and uniform on [0, 2n]. The probability density of x, y, z is therefore 

f(x,y,Z) = ^3 j x z (a 2 -x 2 -z 2 y l/2 x 

5(x - xq cos ft)6(y - ^xq sin ft)6(z - Zo cos <p)dxodzodftdtp. (21) 

The integration can be readily carried out by thinking of xq, ft as a pair of plane polar 
coordinates and transforming to the corresponding cartesian coordinates; similarly with 
Zo, ft- The result is f(x, y,z) = 3 KAn^a 3 ), which is just the inverse of the volume of an ellipsoid 
with semi-axes a, ga, a, as required. 
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Since the evolution of position with time is just a translation in ip and <p, it is clear that 
the space density remains uniform on this ellipsoid. 

The above coordinates and velocities are defined in the cluster-centred frame which moves 
around the galaxy. To transform to the galactocentric rest-frame, we use the equations 

X(t) = (D + x) cos(Qf) + y sin(Qf) (22) 

Y(t) = -(D + x)sin(Qf)+ycos(Qf) (23) 

Z(t) = z (24) 

X(t) = (i + yQ.) cos(Qf) + (y - (D + jc)Q) sin(Qf) (25) 

Y(t) = -(i + yQ) sin(Qf) + (y - (D + x)Q) cos(Qf) (26) 

Z(t) = z (27) 

To ensure that the initial distribution is not bound in the usual sense, its density must lie 
below a certain critical density. Roughly speaking this is the density po of the background, i.e. 
the density which causes the galactic potential, and so 

P< ^= P « (28) 
(cf. eq.O)- 



T he simulations were performed with the particle-mesh code Superbox ( Fellhauer et al 
2000) with 1,000,000 particles and meshes with n — 128 grid points per dimension, to achieve 
high resolution. We chose the following grid sizes for the three nested grids which the code uses: 
^system =12 kpc, R out = 2 kpc and R COK = 0.4 kpc. As the resolution of Superbox is given by 
the length of one grid cell, these parameters give a resolution of 6.45 pc (res = 2R cole /(n - 4)) in 
the central region of the particle distribution. The code has a simple leapfrog integration scheme 
with fixed timestep, which in our simulations was DT= 0.1 Myr. To ensure that our conclusions 
do not depend on the parameters of the code, we repeated one simulation with lower resolution 
and lower particle number, and found no significant change in the results. 



3. Results 

Using initial conditions specified in the previous section, we have conducted a parameter survey 
for various choices of the semi-minor axis a and the density p (Tab. 0. The distance to the 
galactic centre and the circular velocity were held at D = 10 kpc and 220 kms~', respectively. 



3.1. Lifetime of the systems 

In every simulation which satisfied eq. (12 8 1 initially, the distribution was quite stable for several 
orbits around the galactic centre. Fig. [2] shows the Lagrangian radii of the simulation where 
a — 50 pc and p/po = 0.1. The orbital period around the galaxy in our simulations was about 
280 Myr. Clearly the distribution did not evolve much for almost 6 revolutions around the galactic 
centre. 



M. Fellhauer & D.C. Heggie: Unbound Stellar System 7 

Table 1. Parameter space of the simulations. Each simulation performed is marked with a bullet 
point. 



a 
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Fig. 2. Lagrangian radii (10, 20, 90%) for the simulation with a — 50 pc and p/po = 0.1. The 
large figure shows the time interval during which the distribution is stable; the small figure shows 
the whole simulation. 



Because this is an unbound configuration there is no conventional definition for its time of 
dissolution. For the purpose of comparing the simulations, we adopt the time, fdiss, when 

r h = l.5r hfi (29) 

where r/,, r/, o denote the current and initial half-mass radii, respectively. With this definition, the 
dissolution time is shown in Fig.|3]as a function of the density, for different values of a. When p 
is negligible, the particles describe independent epicycles with a common guiding centre, and so 
/diss is essentially infinite (in the epicyclic approximation). 

Except at the highest densities the results can be fitted approximately with power laws 

fdiss=A^J T . (30) 



8 
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p[p ] 



Fig. 3. Dissolution time as a function of the initial density of the particle distribution for different 
values of a. Three-pointed stars correspond to a — 50 pc, crosses to a — 500 pc and five-pointed 
stars to a — 1000 pc. The dissolution time is a decreasing function of p, except forp/po > 1.0. 
In this case the distribution is indeed unstable (according to our criterion), but before expanding 
it undergoes a slight initial collapse which leaves behind a small bound core. 

Table 2. Fitting parameters for the power law relating the dissolution time to the initial density 
of the system. 



a [pc] 


A [Gyr] 


7 


50 


0.86 


0.41 


500 


0.65 


0.35 


1000 


0.60 


0.29 



Best-fitting values for the parameters are shown in Tab.|2] 

We have shown in eq. (1 1 31 that the exponent of the exponential growing solution, k, is ap- 
proximately equal to the square root of Bj which itself is proportional to the square-root of the 
density. This approximation is valid for very small initial density. If, however, we solve ea.(ll3> 
numerically over the range of densities corresponding to our simulations, and fit the results by a 
power law, we obtain approximately 

k cc p 0-40 ± 0.<» (31) 

This shows that, as long as the size a of the system is sufficiently small that the epicyclic approx- 
imation holds, the dissolution time should be almost proportional to p -0, . This is confirmed by 
Tab.|2 
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The effect of the exponentially growing term can be clearly seen in Fig. @J which shows 
particle orbits from the simulation with a — 50 pc and p/po = 0.1 up to about tfo ss . Evidently, the 
exponential instability is much stronger in the y-coordinate than in the x- coordinate; this is to 
be expected, by ea. (ll4> . The result is that the expansion of the system takes place mainly in the 
direction along the orbit around the galactic centre. 

While Fig.|4]shows the effect of the exponential instability on individual orbits, Fig.[5]shows 
its influence in the expansion of the system as a whole. At very early times it is hidden in the 
random noise of the simulation, or indeed by the epicyclic motion itself. Its presence is most 
clear at intermediate times; later on, when the self-gravity of the system becomes negligible, the 
expansion is sub-exponential. 

We also investigated whether there is a relation between the size of the system and the time 
of dissolution, but as one can see in Fig.[6]there is no obvious analytical relation. The fact that 
small systems (small a) have longer dissolution times than large systems can be understood by 
noting that the epicyclic approximation becomes less accurate for larger systems. 

3.2. Appearance of the systems 

Though the initial conditions are highly idealised, it is interesting to look briefly at how the 
appearance of such a system evolves as it disperses. Even after the system starts to disperse, 
forming tidal tails along its orbit, we still find a measurable density enhancement at the position 
of the unbound cluster. In Fig. we show the measured surface density in one of the runs at 
t = 1, 5 and 10 Gyr. The surface density distribution resembles that of an object with a large 
core radius. At late times the surface density is almost flat until nearly 100 pc and then it turns 
over to a shallow power law with an index of approximately -1.6. Fitting an exponential gives 
an exponential scale-length of approximately 100 pc. 

The kinematics of the model depend on the frame in which it is viewed, and are complicated 
by its initial rotation. At t — it can be seen, by setting A — in eas.( ll5t -( ll7t . differentiating, and 
setting t — 0, that i = -*■')>/£ y = t;id x and z = -Zov' sin (p. Thus along the x- and y-directions the 
system appears to rotate, though with different angular velocities. In the z-direction the velocity 
dispersion, <x,, follows the dispersion of zo, i.e. the thickness in the z-direction. We also see that, 
at a given distance r x from the x-axis, the dispersion in x is given by <j 2 = (x 2 ) = K ,2 (y 2 )/^ 2 = 
K ,2 r 2 /(2^ 2 ), with a similar result for <r 2 . As the time of disruption (as defined here) is approached, 
the profiles of cr x and cry approach that for cr z , which itself changes little from that at t — 0. The 
evolution of cr x is faster, however, because the orbits diffuse transverse to this line of sight more 
quickly (Fig.|4}. 

In the galactocentric frame at t — one sees from eas.J22>- (!27> that the corresponding results 
are X = Y(Cl - k'1%) and Y = -Q.D + (X - D){-Q. + %k'), while Z = 1 Again there appears to 
be rotation (with different angular velocities) in the X and Y directions, but for X it can be shown 
from eqs.0 and dill that the angular velocity vanishes if p = 0. As with the velocity dispersions 
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Fig. 4. Upper: Orbits (in the co-moving cluster frame, for t < 2 Gyr) of three particles in the 
simulation with a = 50 pc and p/po =0.1. For clarity two of the orbits have been shifted in 
the x-direction by ±100 pc. One clearly sees the growing deviation from simple epicyclic mo- 
tion. Lower: x (short dashed) and y (solid) residuals (i.e. coordinates with the epicyclic solution 
subtracted) for the above three particles as functions of time. The dotted (x) and long dashed (y) 
lines show fitting functions for an exponential growing solution. Again for clarity the curves of 
two particles are shifted by ±25 pc. 



in the rotating frame, we find that all three velocity dispersions become comparable as the time 
to disruption approaches. Beyond the original radius of the system, however, the profiles of cr x 
and cr y increase with projected radius, because of galactic rotation. 
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The system which is shown in Fig.Qhas a total initial mass of about 660 M Q . At t — 10 Gyr 
the central line-of-sight velocity dispersion is about 1 kms -1 and the core radius is approximately 
100 pc. If one followed the arguments of Mateo (1998) to calculate the virial mass (assuming the 
object is bound), using his formula 

Mat « 167 ^r core of os , (32) 
with [i = 8, the result is roughly 10 5 M . 



4. Conclusions and discussion 

We have studied the evolution of a special class of unbound systems orbiting at a constant dis- 
tance from the centre of an isothermal galaxy model. Initially the systems have uniform density 
within an ellipsoidal region. They move according to epicyclic theory (modified for the self- 
gravity of the system), with guiding centre at the centre of the ellipsoid. These systems are exact 
self-consistent equilibrium models within a steady tidal field. 

While these equilibrium models are believed to be new in the context of a tidally limited stel 



lar system, they have properties in common with Freeman's models of a barred galaxy llFreeman 



1966; 



Binnev & Tremaine 



1987). There too the effective potential inside the system is quadratic 



in a uniformly rotating frame. Freeman's models are the two-dimensional limit of a flattened el- 
lipsoid. It would be interesting to explore this link further. 

Our simulations show that, depending on their initial density and radius, the systems we 
have constructed survive with only modest expansion for a surprisingly long time. When our 
results are scaled to a circular orbit of radius 10 kpc with a circular velocity of 220 kms -1 , a 
system of semi-minor axis 50 pc and density 0.2 times that of the background galaxy survives 
for up to several Gyr (i.e. about 20 Galactic orbits). For small systems the dependence on initial 
density, and the direction in which the dispersal takes place, can be understood on the basis of 
the modified epicyclic theory. For large systems the lifetime may be limited by the accuracy 
of epicyclic theory. Even though the systems gradually disperse, there are measurable density 
enhancements even after a Hubble-time. 

Our results have interesting parallels with those of 



Klevna et al. 



(2003), even though the 



initial conditions and galactic field are completely different. They investigated the survival of an 
unbound clump of particles (with an initial Gaussian velocity distribution) in the inner part of 
the dark halo of a dwarf satellite. As long as the dark matter halo is not cusped, they too found 
that the system disperses only on long time-scales (longer than a Hubble-time). The longevity of 
their models relies partly on the relatively cold initial conditions and partly on the near-harmonic 
potential of the core of the host galaxy. Our models are also cold initially (in the sense that the 
velocity dispersion in x— and y- at a given point is zero), and the potential is approximately 
harmonic inside the model (eqs.(|6j— (JHJ). 

Idealised models, such as those discussed in this paper, have the advantage that they are 
suitable for analytical study, which can provide greater understanding than the compilation of 
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results from simulations. An example of this is our result on the dependence of dissolution time 
on density. The drawback of idealised models is that the results cannot be applied directly to real 
systems, and can only be suggestive. On the other hand a study of Fig|2]suggests that the role of 
the special initial conditions may be confined to the phase in which the Lagrangian radii do not 
evolve by a large factor. The subsequent evolution (e.g. as shown in late-time plots in Fig[7J may 
be applicable more widely. 

The systems in nature which most closely correspond to ours are open clusters. Shortly after 
the tidal radius of an evolving cluster has shrunk to zero, we expect that the remnant consists 
of an inhomogeneous spatial distribution of stars with a preponderance of retrograde rotation, in 
nearly circular motion in the disc of the galaxy. Our models are initially homogeneous, with pure 
retrograde rotation. They nevertheless suggest that it would be interesting to perform similar 
simulations with more realistic initial conditions. If these also exhibited the long dissolution 
times we find for some models, this would suggest that it may be possible to detect the remnants 
of dissolved star clusters even now. 

Other subsystems within galaxies (satellite galaxies and globular clusters) differ further. Their 
galactic orbits are not nearly circular, as a rule, and the resulting time-dependent tide presumably 
causes faster dissolution, which is enhanced also by disk shocking. For this reason, we consider 
that the dissolution time we find (when suitable scaled to the orbital period) is an upper limit to 
the lifetime of such systems. 
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Fig. 5. Upper: Difference between the current half-mass radius r/, and the half-mass radius at the 
beginning of the calculations r^o- All values are taken from simulations with a = 50 pc and vary- 
ing density (from left to right: 1.0, 0.5, 0.1, 0.05 and 0.01 times the background density). Lower: 
Ratio of rh and r^o for the same calculations. (The almost horizontal plot is a test calculation in 
which self-gravity was neglected in both the initial conditions and in the run. For this model the 
only evolutionary mechanism is the inexactness of the epicyclic approximation.) The phase of 
exponential growth is clearly visible. If the half-mass radius gets very large then the self-gravity 
of the system is negligible, and one sees the effects of epicyclic drift: the curves change to a 
linear growth. 
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Fig. 6. Dissolution time as a function of the size a of the system. The density is kept constant at 
p/po = 0.1. 




Fig. 7. Surface density profile of the simulation with a = 50 pc and p/po = 0. 1 at times t = 1 Gyr 
(crosses), t = 5 Gyr (five-pointed stars) and t = 10 Gyr (six-pointed stars) measured in concentric 
radial bins centred on the object and looking along the Z-axis. Also shown are the fitting lines 
for the t = 10 Gyr data. Short dashed line is the horizontal fit for the inner part, long dashed line 
the power law fit for the outer part and the solid curve is an exponential fit. 



